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ABSTRACT 



A hypothesis is made that the Galerkin Finite Element 
Method (GFEM) offers a viable option to the traditional 
Finite Difference Method (FDM) for numerical weather predic- 
tion. The shallow water barotropic primitive equations are 
the forecast equations for all experiments. The hypothesis 
is tested by observing simple, analytic, atmospheric wave 
propagation on uniform and variable mesh grids. Second, a 
strongly forced solution simulating small scale nonlinear 
interactions is evaluated for both the GFEM and FDM. 
Finally, a variable, moving grid for a GFEM model is 
compared to a uniform, higher resolution GFEM model for a 
strong vortex in a mean flow. The GFEM shows a better 
propagation for simple atmospheric waves and better predic- 
tion to a forced nonlinear solution than the FDM model. A 
moving variable grid follows an area of strong gradients 
while not generating noise in the transition zone. 
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I. INTRODUCTION 



Proper simulation of atmospheric flow is the main 
objective of numerical weather prediction. The foundation 
of numerical weather prediction is a set of equations 
including the momentum, thermodynamic, moisture and conti- 
nuity equations for modeling the atmosphere. Through 
computer simulation, numerical weather prediction predicts a 
future state based on initial conditions describing the 
atmosphere. A successful forecast model must include small 
scale processes. Transports, conversions and exchanges of 
mass, momentum and energy occurring on the small scale 
represent important features which must be properly simula- 
ted. Feedbacks from the proper representation of these 
small features can sometimes markedly influence the larger 
scale solutions. Representation of the effects of small 
scale processes can be accomplished directly through 
increased spatial resolution. Continued increases in the 
spatial resolution will eventually allow the desired process 
to be resolved properly. However, a doubling of resolution 
generally requires an eight fold increase in computational 
effort. The value of increased spatial resolution must 
ultimately be measured by its contribution to the overall 
forecast. A level of confidence in the forecast must be 
achieved that the process is resolved near that particular 



14 



grid resolution. The grid resolution could be uniformly 
fine in the forecast domain. The computational effort for 
uniform fine mesh models could be far in excess of available 
resources. An alternative to a uniform fine mesh is a 
variable mesh where the fine mesh covers only regions of 
interest or high activity. 

The conventional numerical weather prediction forecast 
scheme, the finite difference method, approximates the 
partial differential equations with a truncated Taylor 
series. It has performed admirably when forecasting the 
larger scales of motion. Technological improvements in 
computing power coupled with better understanding of the 
atmosphere now allow the smaller scales to be forecast. 

An alternate method for numerical weather prediction, 
the GFEM approxi mates the partial differential equations 
while minimizing the error between the actual equations and 
their approximation. This best fit logically leads one to 
the expectation that the GFEM will better model the smaller 
scales than the finite difference schemes. This research 
will demonstrate practical aspects of the GFEM theoretically 
possi bl e . 

In this dissertation, the GFEM will be evaluated to 
determine its potential to model atmospheric flow. Equiva- 
lent GFEM and FDM models will be utilized to compare the two 
methods. Simulations of small scale processes explicitly 
resolved on uniform and variable grids will be investigated. 
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Rigorous experiments will demonstrate both GFEM and FDM 
responses for forcing near the grid length scale. Finally, 
a demonstration will be made of a GFEM moving variable mesh 
which moves with the small scale process or feature, and 
thereby allows the fine mesh to resolve the highly active 
region initially and throughout the forecast period. 
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II. HYPOTHESIS 



Numerical weather prediction has steadily improved in 
the simulation of atmospheric flow. Large scale flows have 
been adequately represented by finite difference models for 
several decades. Gal erki n- type formulations (Cullen, 1974b; 
Hinsman, 1975; Staniforth and Mitchell, 1978; Cullen and 
Hall, 1979; Staniforth and Daley, 1979; MacPherson and 
Aksel, 1980; and Sasaki and Reddy, 1980) have been shown to 
be competitive with finite difference models, but have not 
shown a marked improvement. Comparing the current opera- 
tional models at selected large computer centers, one finds 
that two are Galerkin and three are finite difference 
(Galerkin: National Meteorological Center and Canadian 

Meteorological Center; and finite difference: Fleet 
Numerical Oceanography Center, Air Force Global Weather 
Center and European Center for Medium Range Weather 
Forecasting. However, all centers have ongoing research 
with Galerkin models and there are indications that these 
will give better long range forecasts. The dichotomy arises 
because the Galerkin applications have not vindicated them- 
selves with a marked increase in accuracy, but rather have 
shown equivalent accuracies. Staniforth and Mitchel (1977) 
stated that ultimately the best gl obal /hemi spheri c models 
will be a spectral model. The spectral model is based on 
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the Galerkin procedure and the use of trigonometric func- 
tions is particularly appealing for hemispheric or global 
grids. However, where non-uniform grids are required, the 
GFEM is a more logical choice. 

Marked improvements in regional forecasts will not come 
from better large-scale models, but rather from models that 
can simulate smaller atmospheric features. The National 
Weather Service Limited Fine Mesh Model is an improvement 
over the hemispheric model, partly because it has higher 
resolution and can resolve smaller phenomena. These smaller 
features can also affect the large scale flow. Numerical 
meteorologists have realized this for years and have 
attempted to model these features. 

The GFEM has the potential to increase efficiently the 
spatial resolution for the purpose of simulating accurately 
the small-scale processes. If a variable mesh is to be 
employed, then it should be evaluated by proper simulation 
of a small-scale feature. In previous research (Staniforth 
and Mitchell, 1978), the refined grid was tested by pro- 
pagating synoptic-scale waves into the finer grid and 
demonstrating that the wave could move into the finer grid 
without generation of significant noise. These conditions 
represent a prerequisite. If synoptic-scale waves cannot 
move freely into and through the variable grid, then advec- 
tion interactions will not occur properly in the fine grid. 



18 



However, one must also show that a smaller scale atmospheric 
feature can be properly resolved or developed in the fine 
mesh area. This will be a milestone. The efforts of 
Staniforth and Mitchell (1978) and Staniforth and Daley 
(1979) have stressed the movement of synoptic scale features 
into the fine mesh area but have not demonstrated improved 
resolution of smaller scale atmospheric features in the fine 
mesh area. 

The hypothesis is that the GFEM is a viable option for 
numerical weather prediction when simulating atmospheric 
flow on variable grids. Three separate features of the GFEM 
will be explored. Each feature will establish the creden- 
tials of the GFEM as a viable option. Each feature is 
intimately related to the proper representation of a small- 
scale phenomenon. The cost effectiveness of the GFEM, when 
evaluated in these contexts, will provide a measure of the 
potential contribution of the method. 

A. FEATURES 

The following three specific features will be explored 
to support the hypothesi s. 

1 . Variable Grid 

Investigations of a suitable alternative to the 
finite difference models for a variable grid will be 
performed. Two basic subdivisions are avai 1 abl e-- tri angul ar 
or rectangular. Staniforth and Mitchell (1978) employed the 



19 



variable rectangular grid as shown in Figure 1. The grid 
contained some areas of finer resolution which were not in 
the verification area. Two points associated with having 
variable resolution in peripheral regions arise: 

(a) Unnecessary computational overhead is required; 
and 

(b) Undesirable phase changes occur as the wave 
propagates in the peripheral regions. 

Older (1981) and Woodward (1981) developed s transformation 

procedure to vary smoothly the resolution for a channel 

domain from a coarse to a fine area in a triangular subdivi 

sion. A uniform equilateral subdivision is shown in Figure 

2. The use of triangles allows the increased resolution to 

occur only where desired and not in peripheral regions. 

Two possible choices of triangles include the right and the 

equilateral. Hinsman (1975) utilized equilateral triangles 

while Cullen (1974b) utilized near-equi 1 ateral triangles. 

reported excellent wave propagation. Kelley and Williams 

(1976) utilized right triangles and experienced very noisy 

solutions. Woodward (1981) duplicated Kelley and Williams 

effort with equilateral triangles and found a major 

reduction in the noise. 

The differing geometries and results thus far reported 
mandate a further review of triangular subdivisions. It is 
not obvious which subdivision is most suitable. A distinct 
advantage of the rectangular subdivision is that it allows 



Both 
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F i 9 • 1* Rectangular subdivision of the Northern Hemisphere 
on a polar stereographi c projection (Staniforth 
and Mitchell, 1978). 
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Fig. 2 . Triangular subdivision for a channel in Cartesian 
coordinates (Woodward, 1981). 
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algorithms to be developed which take full advantage of 
vector processors. However, both subdivisions afford the 
luxury of obtaining a variable grid. Therefore, similar 
experiments performed on each type subdivision will point 
out the advantages and disadvantages of each. 

2 . Numerical Simulation of Physical Processes 

The variable grid chosen as the most suitable 
alternative will be tested while simulating a physical 
process to show that a small-scale atmospheric feature can 
be properly portrayed in the fine mesh area. As stated in 
Chapter II, the ability of variable grids to resolve small 
scale atmospheric features has not been rigorously tested. 

The ability to move synoptic scale features into a finer 
mesh is a prerequisite and must be shown. The main crux of 
the problem is what is happening in the fine mesh area. A 
proper scheme must resolve a small scale feature near the 
smallest grid length. Schoenstadt (1980) and Williams 
(1981) indicated that the most responsive schemes were either 
a staggered finite element grid with primitive variables or 
an unstaggered finite element grid with vorti ci ty/di vergence 
formulation. This research will utilize the latter. 

The simulated process will be that of a mass source 
analogous to that found in the upper atmosphere above a 
hurricane or on the leading side of a strong trough. The 
source will appear in the continuity equation. A known wave 



23 



form will be assumed and the required source to support the 
wave form will be analytically derived. The expression for 
the wave form will allow control of the scale of the process 
so that a near-grid scale phenomena can be simulated. 

3. Moving Grid 

The ability to move the fine grid so that it 
remains centered on an atmospheric circulation will also be 
demonstrated. This capability will further enhance the 
applicability of the GFEM for atmospheric simulation. If it 
is shown that small-scale forcings are properly reflected in 
the flow, and the grid can be moved with the forcing 
disturbance, then the GFEM has great potential for 
atmospheric prediction. 
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III. THE DEVELOPMENT OF THE GFEM 



The hypothesis in Chapter II proclaims the GFEM as a 
viable option for simulating atmospheric flow. An under- 
standing of the history and principles of the GFEM will give 
further insight for the expected improved performance as 
compared to the conventional FDM. 

The origin of the GFEM can be traced to the seventeenth 
century work by Leonhard Euler. His work established the 
branch of mathematics known as calculus of variation. He 
recognized tiat there existed a partial differential equa- 
tion (PDE) associated with the minimization of a functional. 
This PDE has been appropriately named the Eul er-Lagrange 
equation. Solving the PDE was equivalent to determining a 

«k 

stationary value of the functional. In the late nineteenth 
century, Lord Rayleigh furthered the variational calculus by 
representing the dependent variable as a mathematical 
expression, typically as a power series. The stationarity 
of the solution allowed determination of the unknown coeffi- 
cients of the power series. His work was generalized in the 
early twentieth century by W. Ritz and the procedure has 
since been referred to as the Rayl ei gh-Ri tz method. A 
limitation in the solution by the Rayl ei gh-Ri tz method is 
the fullness of the matrix to be inverted. Shortly after 
Ritz completed his work, Galerkin developed a procedure 
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using weighted residuals which had a wider application than 
classical variational calculus. The residual is made ortho- 
gonal to a function, called the test function, and thereby is 
minimized. While the Rayl ei gh-Ri tz procedure applies to the 
functional of an associated Eul er-Lagrange equation, the 
Galerkin procedure pertains to any PDE. The methods produce 
identical results when applied to equivalent extremum 
probl ems . 

The Galerkin procedure was unknowingly well established 
by the end of World War II. The rapid technological 
advances made during and after the war, coupled with the 
emergence of the computer, led the aircraft industry to 
develop new numerical procedures for solving stress problems 
in aircraft design. A successful technique was developed 
and mathematicians realized, long after the fact, that the 
Galerkin procedure was being utilized. 

The GFEM is a commonly used subset of the set of 
Galerkin procedures. After subdivision into a set of ele- 
ments, the domain resembles a completed jigsaw puzzle, and 
hence leads to the terminology "finite elements." Figures 1 
and 2 are samples of such subdivisions. The dependent 
variables of the POE are represented as a linear combination 
of known functions, which are usually low order polynomials. 
The same function is employed as the test function. When 
substituted into the PDE, these approximations leave a 
residual. Minimizing the residual completes the procedure. 
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The functions are globally zero except where the dependent 
variables are defined near each individual node. The pro- 
duct of functions remains zero except where the representing 
and test functions are both non-zero. This makes the method 
attractive for computer implementation. It removes the 
matrix "fullness" problem found with the Rayleigh-Ritz 
approach. The minimization process lies at the heart of the 
expected improved performance. Each term of the POE has 
been simultaneously approximated and the error in those 
approximations has been minimized. 
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IV. MODEL DESCRIPTION 



Two different barotropic shal 1 ow- water models will be 
employed to test the hypothesis in Chapter II. Each model 
will have the identical domain, boundary and initial 
conditions. The difference will lay with the partial 
differential equation approximation. In one model, a 
Galerkin approximation to the partial differential equation 
will be used, while the other model will use a finite 
difference approximation. Two versions of the Galerkin 
model will be evaluated. The first version will have trian- 
gular subdivisions and basis functions, while the other will 
have rectangular subdivisions and basis functions. 

A. GALERKIN FINITE ELEMENT MODEL 

The system of equations referred to as the shallow- 
water equations consists of three equations with three 
forecast variables <j>, u and v. The equations are written 

+ U 11 + V 11 + $ (M + ) = 0 (4-1) 

3 1 3 X sy ^ V 3X 3y' 

+ u iiL + v— - fv + 2 ± . o (4-2) 

at ax ay ax 

+ u i* + v 9 v + fu + a± , 0 (4-3) 
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Here <j> is the geopotential height, 



u is the east/west component of the wind, 
v is the north/south component of the wind, and 
f is the Coriolis parameter 

By expanding <j> into a mean ($) and a deviation ( 4> ' ) the 
equations can be written 



here D is the divergence 

K is the kinetic energy (per unit mass), and 
Q is the absolute vorticity. 

The primes will be dropped for the rest of the paper for 
cl ari ty . 

Cullen and Hall (1979) showed that the accuracy of the 
GFEM solution was better for the vorti ci ty-di vergence 
formulation of the sha 1 1 ow- water equations than for an 
increase in resolution with the primitive formulation. 
Williams and Schoenstadt (1980) noted that staggered 
variable formulation of the primitive equations and the 
unstaggered vorti ci ty-d i vergence formulation gave the best 
treatment of geostrophic adjustment for small-scale 
features . 



+ $D + r~( u<j> ' ) + -^-( v ' ) = 0 



(4-4) 




(4-5) 




(4-6) 
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The vorti ci ty-di vergence form of the shal 1 ow- water 



equations also allows the use of a semi -i mpl i ci t time 
scheme. This scheme artificially slows the propagation 
speed of the fastest gravity waves, which allows a much 
larger time step than one could expect for a normal Courant 
Fredri ch-Lewy (CFL) stability criterion. This scheme thus 
offsets some of the extra computational expense required to 
solve the system of equations assembled at each time step. 

The vorti ci ty/di vergence form of the equations is 



relative vorticity. 

The velocity can be written as the sum of the rotational and 
ir rotational components as 



ft + * D + 97 (ud,) + ‘Ty (v ^ ) = 0 



(4-7) 



|f + ^(uQ) + fy ^ vQ ) = 0 



(4-8) 



|D + + A . £( V Q) + -jy( uQ ) = 0 | 

O 3 V 9 U 

Here is the Laplacian operator and c is the (— - jy) 



(4-9) 



V = V, + V 
~x 



where = K X Vip and = Vy 



The equations can be rewritten using 




and 
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resul ting in 



11 + 

3t 



^v 2 x = - £(u*) - ^ (v<p ) 



(4-10) 



9 n 2 , 

Yt 7 * = 



- &<«*> - 



(4-11) 



(4-12) 



frv 2 x + v 2 <j> = t^-CvQ) - -L(uQ) - -i- iH . _L U< 

3t 9X 3 y v ax ax ay ay 

Equation (4-12) can be further manipulated 

V 2 (.ff + *) - jjCCvQ) - !£) - ~( (uQ) + |£) (4-13) 



3X 



3y 



ay 



The domain of integration is a channel with east-west 
periodicity. The boundary condition at a wall is 

V • N = 0 

where N is an outward pointing normal vector. Along the 
northern and southern walls, the v component is equal to 
zero, so that the v equation of motion (4-3) reduces to 

If = ■ fu 

The zonal and meridional components of the wind can be 
written as non-di vergent and irrotational components as 



and 



u = _ il + iX 

ay ax 



v = M + lx 
ax ay 



Then, along the north/south walls where v equals zero, the 
boundary condition is 



|i ♦ . 0 

ax ay 



(4-14) 
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The above condition is imposed by setting 

tp = a constant 

when solving the vorticity equation and 



(4-15) 




(4-16) 



when solving the divergence equation. This is an overspe- 
cification but (4-14) would be difficult to apply. The 
initial conditions which will be presented later are 
specifically selected to satisfy (4-15) and (4-16). 

The semi -i mpl i ci t scheme is implemented by evaluating 
all the terms on the left hand side of the equations as an 
average at time levels (t + At) and (t - At) or with a 
centered time difference as appropriate. All the terms on 
the right hand side are evaluated at time level t. The 
equations become 



7 2 (ff + *(t+At) + «.(t-4t)) = ^r((vq) - If) - jyUuQ) + If) 



The divergence equation (4-19) can be solved for X at 
(t + At)and substituted into the equation (4-17) to yield 



+ $(v 2 x(t + At) + V 2 x(t-At)) = - j"(uQ) - ^y(vQ) 

’ 2 H ■ - &<“<» - 



(4-17) 



(4-18) 



(4-19) 




( 37 (<t)U) + a“7 Uv)) 
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where the overbar denotes an average of t + At and t - At. 

The system of three equations in three unknowns has been 
reduced to two Poisson equations and one Helmholtz equation 
to be solved at each time step. 

The solution procedure involves solving the <j> equation 
(4-20) for a new £ . The divergence equation (4-19) is then 
solved for £ + and by subtraction for |^-. Finally, the 
vorticity equation (4-18) is solved for |^-. There are two 
options as to which variables will be history-carrying: 
either <|>, i|» and x or <j>, u and v. The choice was made to use 
<(», u and v as history-carrying variables. They are updated 
after each time stop following Staniforth and Mitchell 
(1977) by 

(|> ( t+ A t ) = 2 <j> - 4>(t-At) (4-21) 

u ( t+A t ) = 24t(& - £ If) * U(t-At) (4-22) 

♦ v(t-At) (4-23) 

Implementation of the GFEM is accomplished as described 
in Chapter III. The north and south latitudes are input 
parameters, and the north/south distance is subdivided into 
N equal parts, where N is also an input parameter. The 
east/west increment is a function of the north/south incre- 
ment, as explained in the section describing the individual 
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experiments. The Coriolis parameter will be a constant 
equal to the mid-channel value. 

The triangular and rectangular domain subdivisions will 
be demonstrated in the experiment sections. An appropriate 
approximating function for the rectangular subdivision is a 
bilinear function. In either case, the forecast equations 
in Galerkin form are 



-±t7 )V i 'M ((vQ Vi -n V:» v i 



- / ( ^ ((uQ) j v o + ay K j v j ))v i + + W Uv, jV ¥ i 






- ■ M (uQ) jV v i - / ( ^ (vQ) j v j )v i 

I'^Tt x jVj ♦ ijVjlV. - /C^llv QJjVj - If jVjJlV, 



( 4 - 24 ) 



( 4 - 25 ) 



/ 



t— ((uQ) .V . + — .V . ) )V . 

3y w j j 9y j j ; ; i 



( 4 - 26 ) 



Here the integral sign implies an area integral over the 
domain, the j subscript denotes Einstein summation for the 
dependent variables and the i subscript is the ith nodal 
equati on . 
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The strength of theGFEM is that the spatial deriva- 
tives of the dependent variables become derivatives of the 



known approximating function and, therefore, are known 
exactly at each node. When using piecewise continuous 
linear functions, the first derivatives are piecewise 
discontinuous and second derivatives are not defined. 
Therefore, the second derivatives must be handled in a 
different fashion. 



where the <p implies a line integral along the east/west 
or north/south boundaries. The east/west line integrals 
are zero because of peri odi ci ty, but the north/south line 
integrals add additional terms to the equations and are 
satisfied by the boundary conditions. The Laplacian terms 
become 





e 





S 
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/’Vi ■ -M V j 57 V i ■fjy ' 1 j ^ v i + &Vi 



This integration-by-parts procedure allows the use of a 
linear function while approximating second derivatives. The 

final form of the forecast equations is 

f f ~ 3 V 4 3 V 4 ^ 7 3 V 4 3 V • , 

y J 3 X 3 X J 3 y 3 y $( At ) 2 



/ 



( ^ ((vQ) j v j - K jf^J>> v i 



- /V (uQ) j v j * K o ly j))v i + j &> + i*»)j 

'/^? (t ' 4 t)V j Vl + M (U j (t ' At V + v j* t " 4 t) fr i)V i 

'/VjVi 



f;ii< av. av. . 3 *. av. av., . av.,,, 

‘ y ( ’it' 3 a7 J ax’ * at J ay 3 iy 1 1 * ' ,/ ((uQ) j ay j)v i 



-/ 



((vQ)j $J> V, 



■ / ( ft j a7 j fy* + lt J fy j fy’ * + j v j v i ) = 



/ 



- K j fy j))v i 



-/< 



% lCuQ Vj + K j ay j))V i 



(4-27) 



(4-28) 



(4-29) 
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Here the Helmholtz equation (4-27) has utilized 
V 2 X - DIVERGENCE - + fp- 

and 




along the walls. 

The line integral along the north and south boundary 
has been dropped from the vorticity ecuation (4-28), since 
the value of ij> on the boundaries is a constant and therefore 
the value of pjr is zero. The line integral along the north 
and south boundary has been dropped from the divergence 
equation (4-29), because the value of the normal derivative 
along the north/south walls is zero, es stipulated by the 
boundary conditions. The initial conditions will also 
satisfy the condition that the normal derivative of x is 
zero along the north/south walls. 

B. FINITE DIFFERENCE MODEL 

The comparison model to the GFEM model is an adaptation 
of the staggered, primitive-equation model as described in 
Section 7-4 of Haltiner and Williams (1980). The equations 
are in the flux form and employ Scheme C as shown in Section 
7-3 of Haltiner and Williams (1980). The baroclinic model 
described there has been coded to perform in a barotropic 
mode . 
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Since the GFEM is written in differentiated form, 
vorticity and divergence have been diagnosed from the primi- 
tive variables in the finite difference model to allow 
appropriate comparisons. 

C. NUMERICAL METHODS 

Although no attempt is made to optimize computational 
efficiency in this research model, the GFEM must utilize 
efficient numerical techniques to be considered a viable 
option for numerical weather prediction. Various solution 
procedures are available for the GFEM system of equations. 
Successive over- rel axati on was employed during the early 
research stages. Several problems arose which indicated a 
need for a better solution procedure. The successive over- 
relaxation procedure resulted in a bias if the sweeping 
during each pass was in the same direction. Alternating the 
direction alleviated that particular problem, but the number 
of passes required to achieve a desired level of accuracy 
became prohibitive. A second approach employed a direct 
solver using a Gaussian elimination procedure. The matrices 
from the Galerkin procedure are decomposed into upper and 
lower block tri-diagonal matrices. A preprocessing, repre- 
senting the forward substitution stage, can be done once. 

The back substitution must be done each time a solution is 
desired. The coefficients necessary for this back substitu- 
tion are stored in an efficient manner. The particular 
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algorithm is called a “skyline" solver and is described in 
detail in Bathe and Wilson (1976). Skyline refers to the 
compact method of storing only those coefficients required. 
This method has the desired level of accuracy and a high 
degree of computational efficiency. Staniforth and Mitchell 
(1977) employed another direct solver method called, the 
conjugate- gradient method. This method is very computa- 
tionally efficient, even when the non-constant coefficient 
Helmholtz equation is involved. The theoretical operation 
count indicates that for very large domains with many 
degrees of freedom, the conjugate-gradient method wi'l be 
much better than the other direct solver mentioned ajove. 

The Gauss elimination procedure was utilized instead of the 
conjugate-gradient method for reasons of expediency. 

Numerical integration of the three forecast equations 
involves solving first a Helmholtz equation for d> , second a 
Poisson problem for ip and finally a Poisson problem for x* 
The boundary conditions and the equations are well posed for 
the first two equations. However, the normal derivative of 
x is equal to zero for the last Poisson equation and 
requires special attention. The solution plus a constant is 
also a solution. Since only derivatives of x are required, 
the shape of the field is much more important than the 
actual value of x * To avoid computer round off errors, the 
average value of the terms on the right hand side of (4-29) 
is removed at each time step. This sets the constant to 
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zero and does not allow the solution to grow by a constant 
each iteration. 

The current technology of large mainframe computers 
incorporates vector processing. Utilization of this 
ability can be achieved when the system of equations can be 
vectorized. The GFEM is derived in vector formulation and 
is easily vectorized. 



D. INITIAL CONDITIONS 

1 . Simple Atmospheric Waves 

The initial conditions must allow a relative 
amount of control for the desired input parameters as well 
as satisfy the boundary conditions. The forecast model 
history-carrying variables are <f>, u and v. The analytic 
expression for the streamf uncti on , ip, is 

v = A sin 2 p sin - U(y-y mid ) + p (4-30) 

o 



where A 
W 
L 
n 
0 

y mi d 



amplitude of the perturbati on , 
width of the channel, 
length of the channel, 
wave number, 
mean flow speed, and 
middle point of the channel 
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The first term in the expression is the perturbation. 
The second term represents the north/south slope necessary 
to support a mean flow of 0. The third term is the mean 
depth term. The geopotential height, <j> , is related geostro- 
phically to the streamfuncti on , by 



By use of a trigonometric identity, the streamfuncti on is 
written as 



(4-31) 




(4-32) 



The vorticity is given by 




(4-33) 



where 




2irn 

a 2 L 



The divergence is determined through a linearized form of 
the quasi -geostrophi c divergence equation from Chapter 3-2 
of Haltiner and Williams (1980), in the form 



The right-hand side of this equation becomes 

■f *T f 

cos 32* [(— ^U2ai^Aa2 +-^-Ua 2 3 ^') cos 2a^y - - ^"Ua 2 3 2 '] (4-35) 

Assuming a solution form for divergence on the left-hand 
side 



D = cos °2 x ( C j co s 2a^y + C2) 
the constants and C2 become 

f f f ^ 

c i * - (-^- U2a ] 2A “2 + _ |' Ua 2 3 f ) / ( 4 ai 2 + a 2 2 + ) 



and 



f 3 



0 



2 f. 



Ua 0 4/(a 0 + t” ) 



'2 $ v “2 2 2 $ 



Using v 2 x = D and assuming a solution form for x as 

X = cosa 2 x(C 3 Cos 2 aiy + C 4 ) ( 4 - 36 ) 

then C3 and C4 become 

C3 = - C x /( 4 a^ 2 + a2 2 ) 

and 

C4 = - C 2 / «2 2 
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Observe that ip equals a constant and equals zero along 
the north/south walls. With the expressions for stream- 
function, (4-32), anc velocity potential, (4-36), appropri ate 
formulae for u and v can be derived. 



u = U - sin a£x(Acn sin2a^y + <* 2(03 cos 2a^y + C 4 )) (4-37) 



and 



A 

v = cos a 2 x ( & 22 "( 1 - c 0 s 2ajy) - 2 C 3 a-|Sin 2a-jy) (4-38) 

and <j> is given by equations (4-31) and (4-32). 

2 . Source Ter m 

The addition of a source term into the continuity 
equation will test the resolvability of the GFEM model for a 
source. The source is constructed to represent a small- 
scale meteorological phenomenon to provide a measure of the 
response of the GFEM model to forcing near the smallest grid 
scale. The source initial conditions must be able to 
describe a small-scale feature embedded in a large synoptic 
situation. The initial conditions are derived by choosing a 
desired solution and then back substituting to derive the 
form of the source required. This is a unique approach when 
adding a source term. First, a source, S, is added to the 
continuity equation in the form of 



1 ± + 
at 




li + (3U + 3V.) 

ay vv ax ay' 



S <j> 



(4-39) 
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Following the scale analysis of Chapter 3-2 in Haltiner and 
Williams (1980), a new quasi -geostrophi c potential vorticity 
equation is formed with the source term included. 



9 1 



+ u — + 

Ip 9 X 



% 9y^ ? 



f£ 



) + f S = o 



Solving for S leaves 



S 



( J- _ 11 _a_ + 

' 9 t 9y 9X 




(4-40) 



The technique for testing the source term is to pick a 
streamf uncti on and then solve for the source term which will 
satisfy that streamf uncti on. The particular streamf uncti on 
evaluated will be a steady state solution. The source in 
this case is 






ay 



ay 



The particular form of the streamf uncti on is 



2 - i sin 2 ^) 

* * - U(y-y m1d ) - Asin if e 



w 



(4-41 ) 



(4-42) 



where R is a constant. The use of the exponential allows a 
variation in the scale whereby a small scale feature can be 
embedded in a larger flow pattern. Substitution of this 
particular ip into the expression for the source yields 
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_ 1 / r, , A tt . 2 tt y 

s = f“ Cu T STn w e 



l c • 2 ttx 
R S 1 n L 



\ , 1 ir . • 2 tt x 
Hr L Sin L e 



1 . 2 ttX 

R S,n — 



2 2 2 2 
/ A / TT \ • TT_y / o _ - 2 7T X 1 _ . 2 TT X \ 2 1TV \ 

R L Sin W 3cos L. " R sin L W cos W ^ 



— — sin — — 2 2 

. R L / A / tt \ _ 7ry,47i 2irX 2ir ..._4irX\vx 

+ e ( R ( l> S1n W l T Sln — + RL Sln — >>' 



2 X c • ttX^ 

1 / A ir . / mJiL\ a R L -j- 2tix , 

f - * R L ^f) e sin -7— ) 

0 



.1 Sin 2 

(e R L 5ln <;M 2 cos 2 p 



X c l n 2 tX \ 

R sin L ) 



+ 4A(J) 3 )) (4_43) 

E. NONLINEAR EFFECTS 

Cullen (1974a) discussed the requirement to simulate 
the cascade of energy to unresolved wavelengths. To avoid 
accumulation of energy in the shortest resolvable wave- 
lengths, he developed a two-step method for computing the 
nonlinear terms. Implementation of the two-step method can 
be accomplished efficiently. In the formal development of 
the model equations, there are several nonlinear terms. The 
two-step method involves projecting both variables in a 
nonlinear term into Galerkin space and then performing the 
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multiplication. This procedure insures the best approxima- 
tion to that nonlinear term in a weighted-mean sense. 
Additionally, the computational effort for the two-step 
method is much less than the one-step method which involves 
computation of all the possible interaction coefficients. 

F. STABILITY CRITERION 

The stability criterion for the forecast equations is 



A t < 



1 



M t 

AX 



This is more restrictive than a normal Courant-Fredri ch-Lewy 
CFL) condition. It is made up of the normal CFL criterion 
based on advection plus another effect due to rotation, 
(inertial motion). It provides an upper boundary to the 
maximum allowable time step when using a semi -i mpl i ci t 
scheme. In this research, time steps in excess of an hour 
are easily used. Care must be taken not to exceed the 
stability criterion when larger mean flows are experienced 
or when moving to higher latitudes (larger f ) . 

A more detailed discussion of the derivation of the 
stability criterion can be found in Appendix A. 
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V . 



EXPERIMENTS AND RESULTS 



Experiments will be performed involving the features in 
Chapter II to test the hypothesis. Four separate experi- 
ments will be performed. The first experiment will 
establish the overall performance characteristics of three 
specified models. The next three experiments will address 
each specific feature supporting the hypothesis. 

Standards of comparison must be establishea for each 
experiment. The basic comparison will be achieved through 
harmonic analysis of the initial and forecast fields. The 
harmonic analysis is a double Fourier decomposition. The 
harmonic analysis requires that the data be on a uniform 
grid. Because of the triangular and variable grids, an 
interpolation is necessary. A fifth degree polynomial 
surface fitting routine from the International Mathematics 
and Statistics Library (IMSL) was employed to interpolate 
data on the non-uniform grids to a uniform spacing. First, 
a harmonic analysis is performed along columns of constant x 
to obtain the y structure of the initial conditions. Then a 
harmonic analysis of the amplitudes of each y wave number 
along rows gives a true double harmonic analysis. This 

double transform approach is necessary due to the initial 

2 

conditions. The use of sin --j-p- in the y direction is the 
same as using (1.0 - cos -fjp). The boundary conditions are 
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satisfied by a sine wave in the y direction. Therefore, an 
infinite sine series is needed to represent the cosine 
f uncti on . 

The theoretical movement of each wave is obtained from 
an analysis of the shal 1 ow- water equations. The analysis in 
Section 2-6 of Haltiner and Williams (1980) predicts the 
wave speed as 



C = 



U + 



y 2 + f 2 /gH 



( 5 - 1 ) 



where 

C equal s wave speed , 

H equals height, and 
y x direction wave number. 

By generalizing the streamfuncti on to a two-dimensional wave 
form 



ip = A sin — ^ cos y(x-ct) (5-2) 

and after further manipulation, the wave speed becomes 

C = U(l/(1 + f 0 2 My x 2 + u y 2 )) ) 



where 

$ equals average geopotential height, 
y x equals x direction wave number, and 
u y equals y direction wave number. 



48 



A. EXPERIMENT 1 

Establishing wave propagation characteristics on a 
uniform grid will illuminate the properties of each model. 
Three models will be compared: 

1. A triangular subdivision GFEM model; 

2. A rectangular subdivision GFEM model; and 

3. An equivalent finite difference model. 

The GFEM models will use the differentiated form of the 
shallow-water equations with a semi -i mpl i ci t scheme. The 
finite difference model uses the undifferentiated form on a 
staggered grid with a centered (leapfrog) time scheme. The 
triangular subdivision uses equilateral triangles. The 
relationship between base and height for an equilateral 
triangle is 

BASE = 2 x HEIGHT//! 

The x distance therefore has the same relationship to the y 
di stance 



Ax = 2 Ay//3 

To perform comparisons of similar models, the same x, y 
relationship is retained for the rectangular and finite 
di f f erence model s . 

The basic difference between the rectangular and 
triangular models will be in the approximating polynomials. 
The rectangular polynomials are bilinear while the 
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triangular polynomials are linear. The higher order 
polynomials should provide an increase in accuracy. During 
the integration process, many integrals require evaluation. 
Numerical quadrature could be utilized, however, a more 
efficient method is available through the use of natural 
coordinates. Quadrature with no error can be accomplished 
by formula with this method. A detailed description of the 
natural coordinate method is delineated in Appendix B for 
both triangles and rectangles. 

Diagrams for the grids for the triangular, rectangular 
and finite difference models are shown in Figures 2, 3 and 4 
respectively. The domain is 4896.0 km in the y direction 
and 5653.0 km in the x direction. Each model has 12 incre- 
ments in the x and y directions. This gives 156 degrees of 
freedom. All tests will be for a 48 hour time integration, 
a mean flow of 10 m/s and a mean depth of 1000 meters. A 
perturbation of 1 m/s will be added to the geopotential 
field, which includes the mean height plus the north/south 
slope required for the 10 m/s mean flow. This small pertur- 
bation will focus on the linear aspects of the model's 
capabilities. Experiments described below will evaluate 
larger perturbations, and hence the nonlinear interactions. 
Models are evaluated for wave numbers 1, 2, 3 and 4. 
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Fig. 3. Rectangular uniform subdivision for a channel in 
Cartesian coordinates. 
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Fig. 4. Finite difference uniform grid for a channel in 
Cartesian coordinates. 
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B. RESULTS 1 



The propagation characteri sti cs of the GFEM models and 
the finite difference model on a uniform grid will be 
evaluated. Graphical displays will highlight synoptic 
characteristics. Harmonic analysis will show phase 
propagation and wave amplitude distributions. 

The initial and 48-h forecasts for the conditions 
stated in experiment 1 with wave number one on a triangular 
grid are shown in Figures 5 and 6. The southeast/northwest 
orientation in the forecast fields is due to the sine wave- 
form in the y direction. This waveform with model boundary 
conditions requires an infinite number of y modes. 
Fortunately, the amplitudes of the higher wave numbers are 
small. Consequently, there are only a few important modes. 
The multiple modes cause a skewness in the solutions for all 
three models since the y-modes have different phase speeds. 
Additionally, an analytic expression for the solution 
requires an infinite sum of the modes. For this reason, no 
true solution has been computed to compare with the GFEM and 
FDM model s. 

The divergence equation is the most sensitive equation 
of the three forecast equations. A mean depth of 1000 
meters was chosen because the divergence is large. There- 
fore, a critical evaluation of the sensitivity of the GFEM 
model to this important meteorological parameter can be 
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Fig. 5 
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Initial conditions for the GFEM model with a 
triangular subdivision and wave number one. 
Contour intervals are 600 m^/s^ for geppotential 
height, .2 m/s for u and v, .6x10'° s' * for 
vorticity and .2xl0" 7 s' 1 for divergence. Nodal 
points are denoted by an x. 
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Fig. 6 . 



As in Fig. 
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performed. A divergence field which has propagated 
uniformly with the vorticity field is shown in Figure 6. 
Harmonic analysis (not shown) indicates divergence ampli- 
tudes which initially adjust and then oscillate slightly 
about a mean value of .1270x10"^ s~ ^ . The phases of all 
forecast fields are uniform. Synopti cal 1 y , the forecast 
fields <l>, u, and v indicate a smoothly resolved atmospheric 
wa ve . 

A rectangular subdivision 48-h forecast for wave number 
one is shown in Figure 7. Comparisons between Figures 6 and 
7 show little if any difference which indicates that the 
triangular and rectangular subdivisions give comparable 
results. Harmonic analysis of the divergence fields 
indicates a small oscillation about a mean similar to that 
observed for the triangular subdivision. The 48-h forecast 
from the finite difference model for wave number one is 
shown in Figure 8. The equivalent finite difference model 
also has 12 increments in each direction. Comparisons 
between Figures 6, 7 and 8 show that the finite difference 
model vorticity field lags the GFEM models. Additionally, 
the divergence field has lost its areal extent and developed 
a strong gradient of divergence. 

In each of the wave number one tests there were 12 grid 
points representing the wave. Good propagation should be 
expected from all three of the models for a 12-increment 
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Fig. 7. As in Fig. 6 but a rectangular subdivision 
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As in Fig. 



6 but a FDM model . 
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wave. As the wave number is increased with the same degrees 
of freedom, the wave resolution will be decreased, and 
clearer comparisons can be made. 

Figures 9 and 10 are the initial and 48-h GFEM model 
forecasts for wave number 2, triangular subdivision. 
Forecasts for the rectangular subdivision and the FDM model 
are shown in Figures 11 and 12. The triangular and 
rectangular subdivision forecasts are similar. However, the 
FDM model forecast displays additional divergence and 
vorticity centers near the boundaries. Figures 13 through 
16 are a series similar to Figures 9 through 12 except for 
wave number three. These waves are represented by 4 grid 
points. There are indications in Figure 14 that energy is 
appearing in wave numbers other than wave number 3. 
Comparison of Figures 14 and 15 show that the triangular 
forecast contains more noise than the rectangular forecast. 
Figure 16 shows a divergence field from the FDM model which 
is very different from Figures 14 and 15. There is a 
north/south asymmetry in the u and divergence fields. The u 
cells in the northern regions appear to be expanding while 
the southern cells are decreasing. 

Table lisa composite of the harmonic analysis for 
the triangular, rectangular and finite difference wave 
number 3 case (v component amplitudes). The initial and 
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Fig. 9. As in Fig. 5 except for wave number two. Contour 
intervals are 600 m^/s^ for geopotential height, 
.4 m/sforu and v, .2x10"^ s “ ^ for vorticity and 
.6x10“' s' 1 for divergence. 
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As in Fig. 9 but a 48-h forecast. 
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As in Fig. 



10 but a rectangular 



subdi vision. 



62 




r-flXI5 (METERS I "10 r-flXI5 (METERS) "10 T-flXIS (METERS) 

5. 5 6.S 7.5 8.S 3.5 4.6 5.5 6.6 7.6 8.S 3.S 4.5 S.5 6.6 7.S 



GEOPOTENTIAL HEIGHT 



TRU - 18 
P215221 l 




STREAMLINES 

TRU - 18 




X-RXIS 




VORTICITY 

inu - 18 

F21522U 




3.0 3.0 4.0 

X-RXIS l METERS) 



5.0 

Hltf 



DIVERGENCE 

TRU - 18 
F2 152211 




0.0 l JO a.O 3.0 4.0 5.0 

X-RXIS (METERS) "10* 



Fig. 12. As in Fig. 10 but a FDM model. 
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Fig. 13. 



As in Fig. 5 except for wave number three. 
Contour intervals are 600 m^/s^ for geopotential 
height, .4 m/s for u, 2.0 m/s for v, .6x10-5 s -i 



for vorticity and .2x10"® 




for divergence. 
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As in Fig. 



13 but a 48-h forecast. 
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Fig. 15. As in Fig. 14 but a rectangular subdivision. 
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As in Fig. 14 but a FDM model. 
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48-h v component amplitudes for wave number 3 in the x 
direction and the first six possible y modes are tabulated 
bel ow. 



TABLE 1 



Harmoni c analysis (v 


component 


m/s) - 


wave number 


3 




Triangles 


Rectangl es 


Finite Di f f 


ererce 


Wave 


Initial 


48-H 


Initial 


48-H 


Initial 


43 -H 


1 


2.12 


2.14 


2.20 


2.22 


2.20 


2.20 


2 


.01 


.04 


.01 


.02 


.01 


.09 


3 


.41 


.43 


.44 


.44 


.44 


.42 


4 


.00 


.00 


.00 


.00 


.00 


.00 


5 


.05 


.06 


.06 


.06 


.06 


,05 


6 


.00 


.00 


.00 


.00 


.00 


.00 



The minor differences in the initial state amplitudes 
are due to the differences between the triangular and 
rectangular approximation for the initial conditions. 
However, the 48-h amplitudes show that the rectangular 
version has slightly better preserved the initial state 
amplitudes than either the triangular or finite difference 
models. The finite difference model has transferred more 
amplitude to the other modes, especially in the second y 
mode . 

Figures 17 and 18 are the initial and 48-h 6FEM model 
triangular subdivision forecasts for wave number 4. Wave 
number 4 is a severe test for this model. The energy 
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Fig. 17. As in Fig. 13 except for wave number four. 
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Fig. 18. 



As in Fig. 



17 but a 48-h forecast. 
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transfer to other modes is readily apparent. Additionally, 
both the u component of the wind and the divergence show 
considerable noise at the boundaries. Figure 19 is the 48-h 
forecast for the rectangular model. No energy transfer nor 
noise at the boundaries is apparent in Figure 19. Figure 20 
is the 48-h forecasts for the finite difference model. The 
asymmetry in the wave number three finite difference case 
is magnified for wave number four. The divergence field is 
heavily asymmetric. The u cells have developed a dipole 
from each original center. Comparisons of Figures 18, 19 
and 20 show the superiority of the rectangular forecast. 

While the triangular version is noisy, it still has 
maintained the main features better than the finite 
di fference model . 

Table 2 is similar to Table 1 except for the wave number 
4 case. 

TABLE 2 

Harmonic analysis (v component m/s) - wave number 4 

Triangles Rectangles Finite Difference 



Wave 


Initial 


48-H 


Initial 


48-H 


Initial 


48-H 


1 


2.68 


2.60 


2.94 


2.96 


2.94 


2.96 


2 


.00 


.03 


.01 


.10 


.01 


.05 


3 


.51 


.59 


.58 


.56 


.58 


.60 


4 


.00 


.00 


.00 


.04 


.00 


.05 


5 


.06 


.08 


.08 


.08 


.08 


.08 


6 


.00 


.00 


.00 


.00 


.00 


.01 
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As in Fig. 



18 but a 



rectangular subdivision. 
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Fig. 20. As in Fig. 18 but a FDM model 



73 





The harmonic analysis shows that the triangular version 
transfers more energy in the fifth y mode then does either 
other model. However, the rectangular model transfers more 
energy in the second y mode than does either other model. 
Table 3 compares the divergence for the wave number 4 case. 

TABLE 3 

Harmonic analysis (divergence s"*) - wave number 4 
(All numbers are scaled 10 to the minus 10.) 

Triangles Rectangles Finite Difference 



Wave 


Initial 


48-H 


Initial 


48-H 


Initial 


48-H 


1 


2682.0 


1741 .0 


2929.0 


2179.0 


2801.0 


4462.0 


2 


.0 


140.8 


.0 


207.5 


447.4 


5934.0 


3 


521.3 


391.1 


589.8 


440.7 


545.5 


2314.0 


4 


.0 


20.5 


.0 


81.2 


.0 


923.5 


5 


62.6 


60.5 


84.4 


56.3 


68.7 


1382.0 


6 


.0 


10.5 


.0 


8.4 


.0 


882.0 



The finite difference model has greatly amplified the 
initial divergence. A choice between the triangular or 
rectangular grids version based only on these results would 
be purely subjective. However, the response of both grids 
at this high wave number in terms of divergence is highly 
superior to the FDM. 
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Figure 21 is the phase propagation diagram for the 
three models obtained from the double harmonic analysis. In 
all cases, the percentage propagation is for the first y 
mode and the appropriate x mode. The phase propagation 
diagram indicates that the triangular and rectangular models 
are comparable for the low wave numbers and diverge slightly 
at higher wave numbers. In all cases the GFEM models have 
better wave propagation than the FDM model. 

The analysis performed on the uniform grids has shown 
that each GFEM model performs in a highly satisfactory 
fashion. There are minor differences in the phase 
propagation, while the rectangular version has an advantage 
in the control of energy transfers which manifest themselves 
in the divergence. The finite difference model performs as 
well as the GFEM models for the long waves. However, when 
forecasting the shorter waves, the finite difference model 
does not perform as well as either GFEM model. 

C. EXPERIMENT 2 

This experiment will investigate two options for a 
variable grid. Atmospheric wave propagation on variable 
triangular and rectangular grids will be the test vehicle. 
Since this FDM model has no capability on a variable grid, 
it will not be evaluated during this experiment. The 
comparisons will be between a rectangular and triangular 
GFEM model . 
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Fig. 21. Phase propagation diagram for the triangular, 
rectangular and finite difference models from 
Experiment 1. 
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There are many criteria for determination of a suitable 
variable grid. Both rectangles and triangles allow a 
successful implementation of a variable grid. One 
criterion is the accuracy of resolving the atmospheric waves 
that move through the variable grid. A second criterion is 
the increase i ri resolution. Another consideration is the 
ease with which the grid can be refined. 

Cullen (1979) stressed the importance of a smooth 
variation in element size when increasing the resolution 
from coarse to fine. Older (1981) developed a technique for 
transforming a uniform triangular grid into a variable grid 
when working with a GFEM model. He demonstrated that a 
smoothly varying grid allowed atmospheric waves to propagate 
through the transition zone with much less noise generation 
than an abruptly varying grid. Some finite difference 
models with nested grids, for example Harrison (1973), have 
abruptly varying resolution. A generalization of Older's 
techniques allows a similar transformation for a rectangular 
grid. 

The test cases for experiment 2 will be for a simple 
atmospheric wave. A mean depth of 5 km, a mean flow of 
10 m/s, wave number one and a perturbation of 1.0 m/s are 
empl oyed . 
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D. RESULTS 2 

The ability to refine a grid should be easily achieved. 
The method of transforming from a uniform grid to smoothly 
varying should allow one to choose not only what degree of 
refinement but also where the refinement occurs. Older 
(1981) developed a transformation method based on 

X = x + A cos Kx 



where 

X = the transformed grid 
x = the original grid 
A = a constant 
k = 2ir/L. 

Woodward (1981) modified the transformation to include a 

* . 2 kx 

sin — 

for the longitudinal stretching. A trignometric identity 
is empl oyed yielding 

X = x + A(1 - cos kx) 

This research added the capability to control the location 

of the refinement through 

X = x + A( 1 - cos( kx + <5 ) ) 
a X 

The map factor is defined by 

o X 

~ = 1 + kA sinCkx + s) 
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The maximum and minimum values are 



f£ = 1 + kA , = 1 - kA 

3x max 3x min 



The ratio R of maximum map factor to minimum map factor is 



R 



1+kA 

1-kA 



Solving for A yields 

. R-l 

A " k(R+l ) 

For purposes of this experiment R was chosen to vary from 
1.0 to 4.0. A value of R ~ 4 implies that the minimum X is 
one fourth of the maximum X. The Y transformation is 
performed in a similar fashion except that Y = y + B sin Ly 
i s empl oyed where 

Y = the transformed grid, 
y = the original grid, 

B = a constant, and 
L = 2 it / W 

Placement of the high resolution is accomplished by 
determining the value of 6 required to place the minimum map 
factor as desired. For instance, if the refined grid is 
desired in the middle of the channel (L/2) then <5 would be 
equal to tt/2. 
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Figure 22 is the 48-h forecast with the triangular 
subdivision for a uniform grid. This forecast will act as 
the control for the triangular cases. Figure 23 is the 
triangular subdivision, 48-h forecast with R = 2.0. Figure 
24 is the triangular subdivision, 48-h forecast with 
R = 3.0. Figure 25 is the triangular subdivision, 48-h 
forecast with R = 4.0. All of the 48-h u and v fields 
appear identical. This indicates that there is no degra- 
dation in forecast skill when using the same number of 
degrees of freedom and locating many of those unknowns into 
a refined area. There is no noise apparent in the 
transition regime. However, there are definite differences 
in the divergence field, although the magnitude of the 
maximum/ minimum value of divergence is relatively small 
(l.OxlO -8 s' 1 ). 

Table 4 shows the harmonic analysis amplitudes of the 
geopotential fields for the different triangular cases. The 
relationship between initial and final amplitudes remains 
the same regardless of the degree of variability. The phase 
propagation of the v field in the control case is 97.7%, 
whereas it is 96.7%, 95.4% and 94.1% for R values of 2, 3 
and 4. The variation from the control of the phase propaga- 
tion with grid refinement is under four percent. 
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Fig. 22. As in Fig. 6 for R = 1.0. Contour intervals are 
600 m 2 /s z f or geo potenti al height, .2 m/s for u 
and v, .6x10“ s’ 1 for vorticity and .6xl0“°s“^ 

for di vergence . 
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Fig. 23. 



As in Fig. 22 for R 



2 . 0 . 
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Fig. 24. As in Fig. 22 for R = 3.0. 



83 



Y-flXIS I METERS I hIQT r-ftXIS (METERSI *10f Y-ftXIS ittETCRSI nlfl* 

3.5 4.5 5.5 5.5 7.5 3.5 4.5 S.5 6.5 7.5 3.5 4.S 5.5 8.5 7.5 



GEOPOTENTIAL HEIGHT 



TflU • 48 
T2122115 



STREAMLINES 



THU - 
12122115 




U 

THU - 4B 




TO - 48 
12122115 




VORTICITY 

TO - 48 




DIVERGENCE 

TO - 48 
12122115 




Fig. 25. 



As in Fig. 22 for R 



4.0. 



84 



TABLE 4 



Harmonic Analysis (Geopotential m^/s^) 
Tri ang 1 es 





R = 


1.0 


R = 


2.0 


R = 


3.0 


R = 


4.0 


Wave 


Initial 


48-H 


Initial 


48-H 


Initial 


48-H 


Initial 


48-H 


1 


67.98 


67.16 


67.91 


64.65 


67.91 


64.95 


67.91 


63.70 


2 


.00 


.55 


.00 


.69 


.00 


.99 


.00 


.57 


3 


13.55 


12.90 


13.45 


12.55 


13.41 


13.24 


13.29 


11.16 


4 


.00 


.21 


.00 


.37 


.00 


.54 


.00 


.28 


5 


1 .90 


1 .70 


1.67 


1.30 


1.33 


1.62 


1.32 


1.38 


6 


.00 


.16 


.00 


.23 


.00 


.36 


.00 


.21 



Figure 26 is the control case for the rectangular 
version, 48-h forecast. Figure 27 is the rectangular 
subdivision, 48-h forecast with R = 2.0. Figure 28 is the 
rectangular subdivision, 48-h forecast with R = 3.0. 

Figure 29 is the rectangular subdivision, 48-h forecast 

with R = 4.0. Close inspection reveals no noise in either 
the <j>, u or v fields. No degradation in forecast skill has 
been observed due to increased resolution with a rectangular 
subdivision. All <j>, u and v fields are very similar in 
structure. The main difference in the forecasts lies in the 
divergence field where the magnitudes are small ( 10"^s"^) 
because of the specified mean depth. 
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As in Fig. 7 for R = 1.0. Contour intervals are 
600 m 2 /s 2 for geopotential height, .2 m/s for u 
and v, .6x10"° s"* for vorticity and .6x10"° s~^ 
for divergence. 
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Fig. 29. As in Fig. 26 for R = 4.0. 
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Table 5 is similar to Table 4 except for the 
rectangular cases. Phase propagation for the control case 
v field is 98.0%, whereas it is 98.3%, 98.1%, and 97.9% for 
R values of 2, 3, and 4. The variation from the control 
case of the phase propagation with grid refinement is less 
than half a percent. 

TABLE 5 

Harmonic Analysis (Geopotential m^/s^) 

Rectangl es 





R = 


1.0 


R = 


2.0 


R = 


3.0 


R = 


4.0 


Wave 


Initial 


48-H 


Initial 


48-H 


Initial 


48-H 


Initial 


48-H 


1 


68.03 


67.48 


67.98 


65.26 


68.22 


66.45 


68.14 


66.34 


2 


.01 


1.00 


.03 


.36 


.00 


.85 


.00 


1.22 


3 


13.59 


12.95 


13.39 


13.26 


13.49 


13.81 


13.36 


12.39 


4 


.01 


.46 


.05 


.22 


.01 


.48 


.00 


.43 


5 


1.92 


1.81 


1.57 


1 .64 


1.54 


1.59 


1 .46 


1.19 


6 


.02 


.23 


.07 


.04 


.01 


.30 


.00 


.22 



The tests thus far indicate a successful grid refine- 
ment capability for both triangles and rectangles. There 
are neither drastic impairments nor improvements in either 
phase propagation or amplitude changes for either technique. 
Neither version stands out above the other although the 
rectangular version has less phase propagation variation and 
better amplitude conservation. Both versions allow compara- 
ble grid refinement through the use of the procedure 
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previously described. Both allow an easy and straightfor- 
ward implementation of the procedure. 

However, there is a more fundamental difference between 
the versions. Operators, such as a simple derivative, have 
weighting coefficients associated with finite element appli- 
cations as well as with finite difference schemes. For 
instance, a simple-centered, one-dimensional finite differ- 
ence derivative has weights of 1/2 and -1/2. The finite 
element weights are identical to the finite difference 
weights in this one-dimensional case. The weights can 
similarly be computed for two-dimensional operators such as 
a Laplacian. The quadrature formulas used to determine the 
weights are explained in more depth in Appendix 15. The 
geometry of the triangles affects the weights. For 
instance, the two-dimensional Laplacian applied to a finite 
element application on an equilateral triangular subdivision 
gives weights of 1/6 to all surrounding points, and -1 to 
the center point, as indicated in Figure 30. However, the 
weights change if a triangle with a height equal to one half 
the base is employed as shown in Figure 31. The fact that 
some weights are zero is a strong cause for alarm. If the 
triangles are flattened further, the farthest end point 
weights become negative. This is strictly a function of the 
geometry of the triangle and is a clear sign that use of 
the triangles warrants extreme caution. In fact, Williams 
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Fig. 31. Weights associated with a Laplacian operator for a 
triangular subdivision using triangles with a base 
equal to twice the height. 
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(unpublished notes) has determined that the geometry of the 
triangles dictates that unless a triangle with at least two 
equal sides (isosceles) is employed the boundary conditions 
will not be exactly satisfied. 

E. EXPERIMENT 3 

The capability of the GFEM model to resolve atmospheric 
phenomena near the scale of the smallest grid length will be 
demonstrated. This experiment is in response to Task II of 
the hy pothesi s. 

The ultimate motivation for any "improvement" in a 
model is a better forecast. During the past decade, 
improvements have come, to a certain degree, by increased 
resolution. While it is intuitively obvious that increased 
resolution will improve a forecast, demonstrations have not 
been made to show how the model resolves features near the 
smallest grid length. Staniforth and Mitchell (1973) 
increased the resolution with a variable grid but resolved a 
synoptic scale feature in a 37x37 uniform fine mesh area. 

In this experiment, a source term simulating a mass 
source is added to the continuity equation. An expression 
for the source term in terms of a wave form is derived by 
theoretical means. The wave form allows a very small scale 
feature to be simulated. A detailed description of this 
procedure can be found in Chapter IV-D-2. By assuming the 
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wave is in equilibrium and steady state, the solution and 
source are known for all time when the Rossby number is 
smal 1 . 

The development of the source term expression followed 
traditional quasi -geostrophi c theory. One requirement for 
quasi -geostrophi c theory is that the Rossby number be small. 
In the linear case, the Rossby number will be approximately 
0.1. However, in the nonlinear case the Rossby number is 
approximately 0.4 and the quasi -geostrophi c theory 
assumption is violated. 

The first, test will be on a uniform grid for the 
rectangular GFEM aid finite difference models and the second 
test will be on a variable rectangular grid GFEM model. 

High resolution versions of the uniform test will be run and 
act as the control . A mean depth of 5 km was selected. 
Perturbations of 2.5 m/s and 25 m/s upon a mean flow of 
10 m/s will be evaluated for the uniform grid. This will 
allow both the linear and nonlinear aspects of the different 
models to be observed. Here linearity is implied by the 
smallness of the perturbation amplitude. The source term 
solution is nonlinear and nonlinear interactions do occur 
when the amplitude is small. Each test will be integrated 
for 96 h of forecast time. A perturbation of 25 m/s upon a 
mean flow of 10 m/s will also be evaluated for the variable 
grid. 
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The main thrust of this experiment is to determine the 
grid-scale sensitivity for the GFEM model and the equivalent 
finite difference model. Sensitivity will be measured by 
insertion of a source term into the continuity equation. 
However, the theoretical development presented in Chapter 
IV-D-2 is for the rotational component of the wind. Because 
the initial conditions for u and v are non-di vergent , the 
divergence in the model must increase from a zero initial 
state. This increase will require an adjustment process. 
This serendipity effect will be exploited during the 
compa ri sons . 

The linear cases examined represent a small-scale short 
wave embedded in a long wave pattern. Using a small per- 
turbation for the first tests will establish a level of 
confidence that each model responds to the source term 
creditably. The nonlinear cases represent an active vortex 
in a mean flow. The maximum u component is 35 m/s which is 
hurri cane/ty phoon strength. Physically, the source term in 
the divergence field opposes the advection of the vorticity 
field, and “anchors" the steady-state vorticity solution. 

In either the linear or nonlinear tests, the source 
term is nonlinear. By constraining the perturbation to be 
small, the source term plays a linear role. However, when 
the perturbation is large, the source term allows full 
nonlinear interactions. 



95 



F. RESULTS 3 

1 . Linear Case - Uniform Grids 

Figures 32 and 33 are the initial and 96-h fore- 
casts for the rectangular GFEM model. There are 156 degrees 
of freedom in these forecasts. The initial v field has a 
maximum value of _+2.27 m/s across a separation of 4 incre- 
ments. The initial u field maximum and minimum are 12.5 m/s 
and 7.5 m/s. The source tern is shown in the initial fields 
and remains constant throughout the integration. The diver- 
gence value increased to a steady state solution during the 
first 12 hours with an oscillation that died out after hour 
24. Figure 34 is the equivalent finite difference 96-h 
forecast. Figure 34 demonstrates that the finite difference 
model performs well for these initial conditions. Figure 35 
is the high-resolution 96-h forecast for the rectangular 
GFEM model. Close inspection of Figures 33 and 35 shows 
little difference. The low resolution forecast with 156 
degrees of freedom has converged to the high resolution 
forecast with 600 degrees of freedom. Figure 36 is a graph 
of v component amplitudes at hour 96 as a function ofx and 
y wave numbers for the low and high resolution and finite 
difference models. The high resolution (S253115) and low 
resolution (S153115) forecasts are in excellent agreement. 
However, the FDM (F153115) forecast departs from the control 
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Fig. 32. Initial conditions for the GFEM model with a 

rectangular subdivision and source term added to 
the continuity equation. Resolution is low and 
uniform for a 2.5 m/s pert urbat i on. Contour 
intervals are 600 m^/s^ for geopotential height, 
.5 m/s for u and v, .2xl0" 5 s' 1 for vorticity and 
.2x10" s" 1 for source. 
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Fig. 33. As in Fig. 32 for a 96-h forecast. 
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Fig. 34. 



As in Fig. 33 for a FDM model. 
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35. As in Fig. 33 for a high resolution GFEM model. 
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Fig. 36. V component amplitude as a function of x and y 

wave number at hour 96 for the high ($253115) and 
low (S153115) resolution GFEM models and the 
finite difference (F153115) model and 
perturbation = 2.5 m/s. Contour interval is 
.0 5 m / s. 



101 



even with this relatively linear case. For uniformity of 
comparison, similar ranges of x, y wave numbers will be 
shown in Figures 36 and 41. 

2 . Nonlinear Case - Variable Grids 

Figures 37 and 38 are the initial and 96-h 
forecasts for the rectangular GFEM model with low 
resolution. The initial v field has a maximum of +22.7 m/s 
across a separation of 4 increments, but most of the change 
occurs over 2 increments. The initial u field maximum/ 
minimum is 35 m/s and -15 m/s. The source term as shown 
remains constant during the integrations. The adjustment 
process is evident when viewirg 96-h forecasts. The final v 
field has a maximum of 44.3 m/s, and a minimum of -38.0 m/s. 
The final u field has a maximum of 40.2 m/s and a minimum of 
-29.8 m/s. Figure 39 is the equivalent finite difference 
model forecast of Figure 38 and the marked differences in 
the forecast are readily apparent. The maximum final v 
component is 41.2 m/s and the minimum is -47.3 m/s. The 
maximum final u component is 50.6 m/s and the minimum value 
is -38.6 m/s. Figure 40 is the 96-h forecast for the 
rectangular hi gh- resol uti on GFEM model. Comparisons of 
Figures 38 , 39 and 40 show a convergence of the low- 
resolution GFEM version towards the hi gh- resol uti on solution 
with a much poorer showing for the finite difference model. 
Figure 41 is a graph of v component amplitudes at hour 96 as 
a function of x and y wave numbers for the low and high 
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Fig. 37. As in Fig. 32 for a 25 m/s perturbation. Contour 
intervals are 600 mVs^ for geopotential height, 

5 m/s for u and v, .2x10" s~* for vorticity and 

•2X10" 5 s" 1 for the source. 
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Fig. 38. As in Fig. 37 for a 96-h forecast. 
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Fig. 39 . 



As in Fig. 38 for a FDM model. 
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40. As in Fig. 38 for a high resolution GFEM model. 
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Fig. 41. V component amplitude as a function of x and y 

wave number at hour 96 for the high (S258115) and 
low (R158115) resolution GFEM models and the 
finite difference (F158115) model. Perturba- 
tion = 25.0 m/s. Contour interval is 0.5 m/s. 
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resolution and finite difference model. As in the linear 
case, there is strong agreement between the low and high 
resolution forecasts. The departure of the FDM from the 
control is readily apparent in this highly active case. The 
inability of the FDM to properly handle the nonlinear 
interactions manifest itself as large spurious amplitudes at 
high x and y wave numbers. 

The relatively poor showing of the finite difference 
model for the forced cases shown so far can be improved, to 
a certain degree, by increasing the resolution. However, 
increased resolution requires more computational effort. 
Figure 42 is the 96-h forecast with a perturbation of 25.0 
m/ s for the FDM model. The initial conditions are identical 
tc the nonlinear cases shown in Figure 37. The difference 
in the forecasts is the resolution where there are now 576 
degrees of freedom (24X24). The 96-h forecast is certainly 
better than the 12X12 forecast. However, there is some high 
frequency noise in the forecast. Figure 43 is identical to 
Figure 42 except that there are 1296 degrees of freedom 
(36X36). The high frequency noise is readily apparent and 
is a manifestation of nonlinear aliasing. For this forced 
case the finite difference model is unable to achieve the 
same forecast as the low resolution GFEM model. This result 
will be amplified when a model with the same number of 
degrees of freedom as the low resolution model is employed 
but with a variable grid. 
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Fig. 42. As in Fig. 39 with 576 degrees of freedom. 
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This comparison of the finite element and finite 
difference methods for the source term dramatically high- 
lights the potential improvement when utilizing a Galerkin 
approach. Both graphical displays and harmonic analysis 
leave little doubt as to the superiority of the GFEM model 
over the equivalent finite difference model. 

Adding an extra margin to the already established 
superiority can be accomplished by using the 1 ow- resol uti on 
model with a variable grid. Figure 44 is similar to the 
tests in Figure 38 except that a variable resolution grid 
(R = 2.0) was utilized. A 96-h forecast for a variable 
resolution grid (R = 3.0) is shown in Figure 45. Comparison 
of Figures 38, 40, 44, and 45 show tie convergence of the 
low-resolution uniform solution towards the high resolution 
solution as the degree of variability is increased. Figure 
46 is a graph of geopotential amplitudes at hour 96 versus y 
mode wave number (x wave number one) for the uniform, 
variable, and high resolution models. The R = 2.0 case 
better represents the control than the uniform low 
resolution forecast especially at low wave numbers. 

However, at wave number three and above all models have 
similar differences from the control. 

The harmonic analysis demonstrates an improvement with 
increasing variability but is not conclusive. Another 
method to evaluate the improvement with increased varia- 
bility is to look at the difference charts between the 
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As in Fig. 38 for R = 2.0. 
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Fig. 46. Geopotential amplitude versus x wave number 1 and 
y waves 1-6 at hour 96 for the uniform, variable 
and high resolution GFEM models. Perturbation = 
25.0 m/ s . 
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control and each level of variability. The difference chart 
in Figure 47 is the control 96-h forecast (high resolution) 
minus the low-resolution 95-h forecast. Figure 47 shows a 
maximum difference in the center of the vortex. Figure 48 
is the comparison of the high resolution 96-h forecast and 
the R = 2.0 variable grid $6-h forecast. The maximum 
difference has decreased from the uniform case. Figure 49 
is the comparison of the high resolution 96-h forecast and 
the R = 3.0 variable grid 96-h forecast. The difference 
centers at the right-hand edge in Figures 47, 48 and 49 are 
a result of post processing interpolation errors. Although 
the maximum difference in the center of the vortex is 
greater in Figure 49 than Figure 48, the latitudinal 
difference gradient is less. Root mean square differences 
from the high resolution control are 652.2 m^/s^, 608.8 
m^/s^ and 469.7 m^/s^ for the low resolution, R = 2.0 and 
R = 3.0 tests respectively. The forecasts have definitely 
been improved at each level of increased resolution. 

G. EXPERIMENT 4 

This experiment will demonstrate a moving grid 
capability. Harrison (1973) demonstrated the use of a 
moving grid in his finite difference model. The inherent 
problems with moving finite difference models are abrupt 
changes in resolution causing spurious noise. The noise, 
unless damped, is magnified if left unattended every time a 
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Fig. 47. Forecasts and difference charts at hour 96 for the 
uniform high (S258115) and low ( SI 58 1 1 5 ) resolu- 
tion GFEM models. Perturbation = 25.0 m/s. 

Contour i nterval i s 600 m^/s^ for the geopotential 
fields and 200 m /s for the difference field. 
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As in Fig. 47 for the high 
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Fig. 49. 
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As in Fig. 47 for the high resolution (S258115) 
and the variable (S198115) (R = 3.0) GFEM models. 
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grid moves. In a GFEM model there is a smooth transition 
into the fine grid area with no attendant noise. Therefore, 
moving the entire grid with the fine grid is simply a matter 
of accurate interpolation. 

Three tests will be performed. First, uniform low and 
high resolution models will be integrated. These models 
will not use a moving grid. Then a model with the same 
degrees of freedom as the 1 ow- resol Jt i on model, but 
employing a moving, variable grid, will be run. The high- 
resolution model will act as the control. Identical initial 
conditions will be used for all three tests. These condi- 
tions are similar to those developed for experiment 3, 
except that the source term will be zero throughout the 
course of the integration. The initial conditions include a 
sharp trough with a relative vorticity maximum which is 
advected by the mean flow in the absence of any forcing (the 
source) . 

Of particular concern is the generation of noise 
incurred during grid movement. Therefore, a forecast of 
96-h is made during which this grid moved eleven times. As 
previously stated, the model utilizes absolutely no diffu- 
sion, filters or friction. If spurious modes are excited 
during the movement of the grid, they will appear in the 
harmonic analysis. 
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The procedure of moving the finite element grid is 
achieved by displacing all the nodes one unit distance in 
the direction desired. This retains the relative 
distribution of the nodes and merely moves the entire grid 
as an entity. Ramifications of this approach are: 

1. The areas of each element remain the same and the 
numerical quadrature coefficients for all the 
associated integrals do not have to be recomputed. 

2. The unit distance moves should be the smallest 
grid increment; and 

3. Direct solvers that utilize a preproces si ng 
procedure based on the grid distribution will not 
require the preprocessing to be recomputed after 
each move. 

The most important aspect of moving the grid is to accurate- 
ly interpolate the new field values for the new grid point 
position. In this dissertation a standard International 
Mathematics and Statistical Library's (IMSL) bi-cubic spline 
interpolating routine was employed. The movement will occur 
only in the x direction and, therefore, a one-dimensional 
interpolator suffices. The interpolator accommodates the 
periodic boundary condition. However, interpolators for 
non-periodic boundary conditions are available. The accu- 
racy of the interpolator is exact at grid points for a 
uniform distribution. Initial experimentation verified that 
a field could be moved with no loss in accuracy for a 
uniform grid. Experiment 4 will test its capability to move 
a variable grid with R = 2.0. 
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H. RESULTS 4 



Figures 50 and 51 are the initial and 96-h forecasts 
for the uniform low-resolution model. The perturbation is 
25.0 m/s, mean flow 10 m/s and mean depth 5 km. The initial 
divergence is a zero field but spins up through geostrophic 
adjustment. The northeast/ southwest trough orientation is a 
manifestation of the nonlinear interactions occurr - ng for 
this highly active perturbati on. Figure 52 is the 96-h 
forecast for the uniform hi gh- resol uti on model. The 
northeast/southwest trough orientation is evident in this 
forecast as well. Figure 53 is the 96-h forecast r or the 
low-resolution, moving variable grid (R = 2.0) model. 

Time steps were 4320s, 3085.7s and 3600s for :he low, 
high and variable models, respectively. The variable model 
took 95 full time steps and two half steps during the 96-h 
forecast. The half steps are required only to sta'~t the 
model to obtain the N and N + l time levels in the leapfrog 
time stepping. The grid moved 11 times during the forecast 
at time steps 8, 16 , 24, 32, 40, 48, 56, 64 , 72, 80 and 88. 

The x axis origin in Figure 53 is 3530 km indicating 
eastward grid movement. The vortex is still in the fine 
mesh area. Additionally, there is no noise in the forecast 
fields. To reiterate a previous statement, there are no 
filters, diffusion or friction terms employed during this 
integration. Comparisons among Figures 51, 52 and 53 show 
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Fig. 51. As in Fig. 50 for a 96-h forecast. 
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Fig. 52. As in Fig. 51 for a high resolution GFEM model. 
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Fig. 53. As in Fig. 51 for a moving variable (R = 2.0), 
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that the variable grid forecast more closely approaches the 
control than does the uniform low-resolution model. This is 
to be expected since the fine mesh area has been moved to 
coincide with the active region of the forecast domain. 

Table 6 shows the geopotential harmonic analysis for 
the three models after a 96-h forecast. 



TABLE 6 

Harmonic Analysis; (Geopotential) for Experiment 4 

Units are m^/s^ 



Wave 


Low Resolution 


High Resolution 


Moving Grid 


1 


876.1 


904.7 


850.0 


2 


167.4 


152.9 


146.6 


3 


122.5 


176.9 


159.6 


4 


54.7 


39.0 


39.0 


5 


20.3 


10.2 


9.6 


6 


5.1 


7.6 


7.3 



Comparisons of the higher wave number amplitudes show 
almost identical amplitudes in y modes four, five and six 
for the moving case arid hi gh- resol uti on uniform case. The 
amplitude in y mode six indicates that the movement of the 
grid has not generated unwanted noise. The overall improve- 
ment for the moving case is a result of the finer grid's 
ability to better resolve the small-scale interactions 
occurri ng . 
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VI. SUMMARY AND RECOMMENDATIONS 



Numerical weather prediction requires that the 
equations governing atmospheric flow be solved numerically. 
These flows cover a large range of scales. Classical 
modeling schemes have successfully forecast the longer 
scales and improvements are being sought for the smaller 
scales. 

In this dissertation, the GFEM has been employed for 
the shal 1 ow- water equations in a channel domain. The 
differentiated form of the equations on an unstaggered 
finite element grid with a semi -i mpl i ci t time scheme has 
been utilized. An equivalent finite difference model 
representing the conventional approach has also been 
utilized. Analytic initial conditions representing simple 
atmospheric waves and a more complex source term have been 
used to test the different models. 

The hypothesis of this dissertation is that the GFEM is 
a viable option for numerical weather prediction. The 
theoretical foundation of the GFEM leads to a minimization 
of the error between the actual equations and their approxi- 
mation. This minimization implies a basis for expected 
improvement in forecast capabilities. Three features were 
addressed to support that hypothesis. First, alternatives 
for variable grids were investigated. Second, forcing near 
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the smallest grid scale determined if the GFEM model 
resolved the smaller scales better than a finite difference 
counterpart. Finally, a proof of concept for the ability to 
move the grid was demonstrated. Each task adds to the 
usability of the GFEM and illuminates the well-rounded 
potential that it contains. 

First, the ease with which the method allows variable 
grids was demonstrated. In all cases, no numerical 
techniques were employed to damp or filter any field. The 
grid refinement obtained when using a rectangular 
subdivision was more attractive than with a triangular 
subdivision. It allowed refinement to be easily obtained 
with increased accuracy from the higher order polynomial and 
decreased computational effort due to reduced operation 
counts. Additionally, it was shown that the geometry of the 
triangular elements greatly influences the problem. 
Furthermore unless specific constraints are enforced, the 
boundary conditions would not be satisfied. 

Second, the ability of the GFEM to resolve atmospheric 
phenomena occurring near the scale of the smallest grid- 
length was shown. A small-scale steady state solution was 
cast into a source required for that solution. Integrations 
showed that the GFEM model resolved the small-scale forcing 
admirably and far exceeded the performance of an equivalent 
finite difference model. Even when the resolution of the 
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finite difference model was doubled and tripled, it could 
not match the performance of the 1 ow- resol uti on variable 
grid GFEM model. The responsiveness of the GFEM model near 
the smallest resolvable grid length clearly demonstrates 
that the variable grid concept is worth the computational 
effort . 

Finally, the concept of moving the variable grid was 
demonstrated. The ability to resolve a small-scale 
phenomenon is a tremendous asset to the GFEM. The ability 
to move the grid with the atmospheric phenomena further 
enhances its usefulness. 

Throughout this dissertation the GFEM model has 
demonstrated desirable characteri sti cs. It has conservative 
properties. It propagates atmospheric waves better than an 
equivalent finite difference model. It allows variable- 
resolution grids and responds better than an equivalent 
finite difference model near the smallest gridlength. 

Moving grids can be achieved with no apparent noise 
generation. It can utilize direct solvers and is a natural 
choice for vectori zati on on large computers. It truly is a 
viable option for simulation of atmospheric flows. 

The success of this dissertation mandates further 
research. The next logical step would be a baroclinic 
model. Completion of a baroclinic model should allow its 
application as a regional model. Additionally, different 
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basis functions could be employed for the vertical 
discretization, such as solutions to the vertical structure 
equation. For the moving vari abl e-resol uti on grid, further 
research into two-dimensional movement is required. The 
proof of concept in this dissertation of one-dimensional 
movement can logically be extended to two dimensions. Upon 
completion of two-dimensional movement, an operational 
forecast capability for tropical cyclones should be investi- 
gated. The superior small-scale response of the GFEM 
indicates potential increase in skill for both regional and 
tropical cyclone/hurricane prediction. 
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APPENDIX A 



STABILITY ANALYSIS 



A one-dimensional linear analysis of the forecast 
equations will be presented. The analysis was obtained 
through personal communication with Dr. Andrew Staniforth 
(Recherche en Prevision Numerique). The primitive form of 
the forecast equations and a semi -i mpl i ci t time scheme will 
be analyzed, because the results are identical for the 
vorti ci ty/di vergence form. The one-dimensional equations 
with a mean flow U are 



3U 

at 



3V. 

at 



at 



11 = - u 

ax 


|i+ fv 

3x 


(A-l) 


= - U 


11- fu 

ax. 


( A-2 ) 


, .au _ 

+ — - - 

3X 


u 11 

ax 


(A-3) 



Time derivatives are evaluated with a centered time dif- 
ferencing and the other terms on the left-hand side are 
averaged between time levels (t + At) and (t - At). 
Equations (A-l), ( A- 2 ) and (A-3) become 

u(x,t+At ) - u(x,t-At) 1 r (<j> (x+AX , t+At ) - <f> (x-AX , t+At ) 

2At 2 L 2ax 

+ (<j>(x+AX,t-At) - d>(x-AX,t-At)- i _ _ f. r u(x+AX,t) - u(x-AX,th 
2AX J - - L 2 AX -* 

+ f v(x,t) (A-4) 
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v(x,t+At) - V(x,t-At) _ rV(.X+AX,t) “ v(x-AX,t)n 
2At - - U L 2ax j 



- f u(x,t) 



<j>(x,t+At) - <j>(x,t-At) . $ p(u(x+AX,t+At) - u(x-AX,t+At) 
2At 2 L 2 AX 



(u(x+AX,t-At) - u(x-AX,t-At))n _ r. / <i> U+AX , t ) - <f>(x-AX,th 
2AX J " " U( ‘ 2 AX ’ 



Assuming that a function, F, can be expressed as 
F ( x , t ) = F'e 1 ^* + wt > 



then 



F(x,t+At) - F(x,t-At) _ i r 1 

At 



i (kx + wt) 
sin( wAt)e 



F(x,t+At) + F(x,t-At) 
2 

Equations (A-4), (A-5), and 
substitution. 



i (kx + wt) 

= F' cos(wAt)e 



(A-6) - become, after 



isu 1 

At 

fu' + 
ikc$u‘ 



fv' + i kc4> 
isv 1 

At 

is±' 

At 



0 



0 



(A-5) 
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( A- 1 0 ) 

(A- 1 1 ) 
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where s = sin wAt + k'uAt, c = cos wAt and k' = sin — 

The determinant of the system of equations (A-10), ( A - 1 1) 
and ( A- 1 2) when set equal to zero yields 



it ^2 - l? + 



'? 9 

kcl)] = 0 



( A- 1 3 ) 



The roots are s=0 and 



s 2 = (fAt) 2 + c 2 CkAt) 2 $ (A- 14) 

The roots of this second order equation when requiring w to 
be real yield the stability criterion 

At < — J ( A- 1 5 ) 



133 



APPENDIX B 



NUMERICAL QUADRATURE 



Evaluation of the Galerkin integrals requires a fast, 
efficient method . Two different elements, triangles and 
rectangles, have been employed and the quadrature schemes 
for evaluation are presented here. In the first case, 
triangles, area coordinates are used. For the rectangles, 
integration formulas are based on an orthogonal axis trans 
formation. In both cases, the integrals to be evaluated 
contain either products of basis functions, products of 
derivatives of basis functions, or a mixture of both. 

Zienkiewicz (1971) developed the relationship between 
the Cartesian coordinates of a triangle and the area 
coordinates as 



where (x^yj), (x 2 ,y 2 ) and (x 3 ,y 3 ) are the Cartesian coordinates 
of the triangle's vertices and , L 2 and L 3 are the area 
coordinates. Here 



(B-l) 



x = L i x i + l 2 x 2 + i . 3 x 3 




( B-2 ) 




(B-3) 



Li * A, /A 



L 2 = A 2 /A 



L 3 = V A 
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where A ^ , A 2 and A 3 are shown in Figure 54 and A is the area of 
the entire triangle. He further shows that 





dxdy = 



a!b!c! 

(a+b+c+ 2 }! 



2A 



The L's are the basis and test functions used during the 
approximating process. Figure 55 defines some necessary 
parameters for evaluation of the derivatives. 

Equations B-l, B-2 and B-3 can be written in matrix 
form as 
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Eq. (B-5) 
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shows that 



( B-6 ) 



(B-7) 



However, the derivative of the basis function with respect 
to the area coordinate is non-zero only where the basis 
function and area coordinate coincide. In fact, the non- 
zero value is one. Therefore, using Vj as a basis function 



9V 



1 



9X 




(B- 8 ) 
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Fig. 54. Cartesian coordinates vs. natural coordinates. 
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All integrals can be evaluated by combination of these 
quadrature formulas for triangles. 

An orthogonal axis transformation will allow similar 
quadrature formulas for rectangles. The rectangle shown in 
Figure 56 can be transformed by 



c = 



X-X_ 




(B-9) 



The values of c and n at each corner are shown in 

parentheses. The basis function, , becomes 

tl+C i c)(l+n i n) (b-10) 
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Fig. 56. Orthogonal axis transformation for rectangular 
integration formulas. 



Derivatives of the basis function become 
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By substituting and integrating, the 4 x 4 matrix with the 
interaction coefficients can be determined for a derivative 
or straight inner product. For instance, the straight inner 
product is 

(f 1 1 

c ij =//v jdxdy = Sb J ^ dcdn 
1 1 

T6 O+n^r )(1+Tijn)dn 



? 6 ^ 2+ f c 1 c j ^ ^ 2+ § n 1 ^ 



( B -1 2) 



and the mixed derivative is 

•3V. aV. 



1 1 aV . aV. 



ff d'i. 3V . h r f j V . 3V. 

d ij u Jhri^ dxdy = a l l dcdn 



b 

16a 



I I 

J ^ (l + Ti^n){l + nj.ri)d!i 



' il? (2c 1 c j )(2+ f Vj 1 (B ' 13) 

These quadrature rules allow the evaluation of any integrals 
when using rectangles. 
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